Cluster simulations of loop models on two-dimensional lattices 
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We develop cluster algorithms for a broad class of loop models on two-dimensional lattices, in- 
cluding several standard 0(n) loop models at n > 1. We show that our algorithm has little or no 
critical slowing-down when 1 < n < 2. We use this algorithm to investigate the honeycomb-lattice 
0(n) loop model, for which we determine several new critical exponents, and a square-lattice 0(n) 
loop model, for which we obtain new information on the phase diagram. 
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From the beginning of the theory of critical phenom- 
ena, two models have played a central role: the g-state 
Potts model 0, Q and the O(n) spin model [3, 0] . The 
parameter q or n is initially a positive integer, but the 
Fortuin-Kasteleyn (FK) representation [f| and the loop 
representation [6| show, respectively, how the models can 
be extended to arbitrary real or even complex values of q 
and n Q. In particular, for q, n > the extended model 
has a probabilistic interpretation as a model of random 
geometric objects: clusters [§| or loops respectively. 
These geometric models play a major role in recent devel- 
opments of conformal field theory [lfj] via their connec- 
tion with stochastic Loewner evolution (SLE) [ill H3|. 

Since nontrivial models of statistical mechanics are 
rarely exactly soluble, Monte Carlo simulations have 
been an important tool for obtaining information on 
phase diagrams and critical exponents [13]. Unfortu- 
nately, Monte Carlo simulations typically suffer from se- 
vere critical slowing-down, so that the computational ef- 
ficiency tends rapidly to zero as the critical point is ap- 
proached [141 ]. An important advance was made in 1987 
with the invention of the Swendsen-Wang (SW) algo- 
rithm [151 ] for simulating the ferromagnetic Potts model 
at positive integer q, based on passing back and forth 
between the spin and FK representations. The SW al- 
gorithm does not eliminate critical slowing-down, but it 
radically reduces it [lij . Since then, many similar "clus- 
ter algorithms" have been devised, based on this principle 
[13] of augmenting the original spin model with auxiliary 
variables and then passing back and forth. But cluster 
algorithms have traditionally been limited to integer q, 
since they make essential use of the spin representation. 

This limitation was first overcome in 1998 by Chayes 
and Machta [l8j] , who devised a cluster algorithm for sim- 
ulating the FK random-cluster model at any real q > 1. 
For loop models, by contrast, efficient simulation at non- 
integer n has remained out of reach; to our knowledge 
only two Monte Carlo simulations at n ^ 1 have ever 



been published [13,111, and they used local algorithms. 
(Instead, numerical transfer-matrix techniques have been 
employed (ll|.) As a result, many open questions remain: 
for instance, the nature of the phase transition is un- 
clear for the n > 3/2 honeycomb-lattice loop model with 
vacancies [13]; an d the phase diagrams and universality 
classes of loop models on lattices other than honeycomb 
are largely unexplored 23 1 . 

In this Letter we shall set forth a broad (but non- 
specific) generalization of the Chayes-Machta idea, and 
then show how it can be adapted to provide a cluster al- 
gorithm for simulating loop models on two-dimensional 
lattices at any real n > 1. We shall furthermore present 
numerical evidence that this cluster algorithm has little 
or no critical slowing-down when 1 < n < 2. Finally, 
we shall use this algorithm to obtain new results for the 
phase diagram of a loop model on the square lattice. 

We begin by considering a generalized random- cluster 
(RC) model, defined as follows: Let G — (V, E) be a finite 
graph with vertex set V and edge set E. A configuration 
of the model is specified by a subset A C E (the set of 
"occupied bonds" ) , and the partition function is 



Z = 




tlW(Hi) 



(1) 



ACE \e£A 



where H± , . . . , Hk are the connected components of the 
graph (V, A); here {v e } are nonnegative edge weights, 
and {W(H)} are nonnegative weights associated to the 
connected subgraphs H of G. The model (Q} reduces to 
the FK random-cluster model if W(H) = q for all H; 
other special cases include an FK representation for the 



Potts model in a magnetic field, W(H ) = q — 1 + e 



Ji\V(H)\ 

24(, and the loop models to be discussed below. 

Now let m be a positive integer, and let us decom- 
pose each weight W(H) into m nonnegative pieces, any 
way we like: W(H) = YZ=i w <*{H). The first step of 
our generalized Chayes-Machta algorithm, given a bond 



2 



configuration A, is to choose, independently for each con- 
nected component Hi, a "color" a S {1, ...,m} with 
probabilities W Q (/fj)/W(ifi); this color is then assigned 
to all the vertices of Hi. The vertex set V is thus parti- 
tioned as V — UI!Li ^o* ^ ^ s n0 * nar d to see that, condi- 
tioning on this decomposition, the bond configuration is 
nothing other than a generalized RC model with weights 
{W a (H)} on the induced subgraph G[V Q ], independently 
for each a. 

We now have the right to update these generalized RC 
models by any valid Monte Carlo algorithm. One valid 
update is "do nothing"; this corres pon ds to the "inac- 
tive" colors of Chayes and Machta [18(. Of course, we 
must also include at least one nontrivial update! The ba- 
sic idea is to have at least one color for which the weights 
W a (H) are "easy" to simulate. The original Chayes- 
Machta choice, when W(H) = q for all H, is to take 
W a (H) = 1 for one or more colors a (the so-called "ac- 
tive" colors); the corresponding model on G[V Q ] is then 
independent bond percolation, which can be trivially up- 
dated. Since we must have W a (H) < W(H), this works 
whenever q > 1. 

We hope that this explication/extension of the 
Chayes-Machta framework will inspire other researchers 
to invent diverse algorithms of Chayes-Machta type. 
Here we would like to exhibit two such families of al- 
gorithms, for simulating a general class of "loop models" 
on two-dimensional lattices. More precisely, the models 
we have in mind should be called Eulerian- subgraph mod- 
els, because the connected components are not necessar- 
ily loops; we recall that a bond configuration A is called 
Eulerian if every vertex has even degree (i.e., an even 
number of incident occupied bonds; zero is allowed). So 
a generalized RC model is an Eulerian-subgraph model 
if W(H) = whenever H is not Eulerian. The simplest 
such model (the "standard Eulerian-subgraph model") 
has 



as a loop representation of an n-component spin model, 



W(H) = { 



n if H is Eulerian 
otherwise 



(2) 



Another such model (the "disjoint-loop model") has 
W(H) = { 



n 




if H is a cycle or an isolated vertex 
otherwise 



(3) 

These two models are identical if the underlying graph 
G has maximum degree < 3 (e.g. the honeycomb lat- 
tice) but are different otherwise. More generally, we can 
consider the "degree-weighted model" 



W(H) = n J] t dH(i) 

t£V(H) 



(4) 



where djj(i) is the degree of the vertex i in the graph H, 
and to,ti,t2, ■ ■ • are nonnegative weights satisfying t c i = 
for all odd degrees d. 

It is well known [|| that on any graph G of maximum 
degree < 3, the model (H}/© arises for positive integer n 



Z 



Tr J! (1 



j cr l ■ a-j) 



(5) 



where cr = (a 1 ,..., a 11 ) <E K™ and Tr denotes normal- 
ized integration with respect to any a priori measure 
(•)o on R™ satisfying (cr a a^) = n^S* 13 and (a a ) Q = 
(cr Q cr /3 cr 7 )o = 0. (In particular, uniform measure on the 
unit sphere is allowed, as are various "face-cubic" and 
"corner-cubic" measures Q.) For n^l, the Boltzmann 
weight with spins on a sphere defines a non-standard 
0(n) spin model, which has positive weights only for 
\xij\ < 1/n, but it is nevertheless expected to belong 
to the usual 0(n) universality class. 

Likewise, on any graph G, the model [but not f3|] 
arises from the spin model with a "face-cubic" a priori 
measure, i.e. cr is a unit vector ±e Q (1 < a < n) with 
probability l/2n. 

Finally, let us observe that, on any planar graph G, 
there is a one-to-two correspondence between Eulerian 
bond configurations A on G and Ising configurations {a} 
on the dual graph G* (namely, A is the Peierls contour 
for {cr}). Under this mapping, the bond model ([2]) on G 
with n = 1 is isomorphic to the Ising model on G* with 
couplings J e satisfying v e — e~ 2J ' . More generally, the 
model |4]) with n = 1 and tj, — ab d , i.e. 



W(H) = 



a \V(H)\ b 2\E(H) 





if H is Eulerian 
otherwise 



(6) 



gives an Ising model with b 2 v e 

With this observation, we can present two algorithms 
of Chayes-Machta type for simulating Eulerian-subgraph 
models on a planar graph G. In these algorithms, we 
keep at all stages an Eulerian bond configuration A on 
G together with one of the corresponding Ising spin con- 
figurations {a} on G*. 

The first algorithm applies to those models in which 
there exist a, b > such that W(H) > a \ v ^b^ E ^ for 
all Eulerian H . This includes in particular the model ^ 
with n > 1 and t<i > ab d . 

Algorithm 1: We use two colors: an "active" color 
(a = 1) with W\(H) given by ([6]), and an "inactive" color 
(a = 2) carrying the remaining weight. The first step of 
the algorithm leads to a partitioning V = V\ U V2. We 
now freeze all bonds having one or both vertices in V% in 
their current state, while leaving bonds within Vi free to 
move. The latter bonds are updated by applying one step 
of the Swendsen-Wang algorithm to the correspondingly 
constrained Ising model on G*. In detail, the algorithm 
proceeds as follows (see [25( for a proof of validity): 

1) Independently for each component Hi, color it a = 1 
with probability W\{Hi) /W{Hj), and a — 2 otherwise. 

2) On each edge of G* whose dual does not lie entirely 
within Vi, place a bond. On each edge ij of G* whose 
dual e lies entirely within V± and for which J e aiO~j > 0, 



3 



place a bond with probability 1 — e 2 I ,7 =I. (Here J e is 
defined by b 2 v e = e~ 2Jc .) 

3) Form clusters on G* from sites connected by bonds. 
Independently on each cluster, flip the Ising spins with 
probability 1/2. (Note that a cluster typically contains 
spins of both signs.) 

4) The new bond configuration is the Peierls contour 
for the new Ising configuration. 

This algorithm is easily extended to allowing more 
than one active color, in case W{H) > fca |l/(ff)l 6 2|B(H)l 
for some integer k > 2. We freeze all bonds having one 
or both vertices in an inactive color as well as all bonds 
connecting two distinct active colors. 

Please note that Algorithm 1 is actually a family 
of algorithms corresponding to different choices of a, b. 
Among the maximal allowed pairs (a, 6), some choices 
may be more efficient than others. 

Alternatively, one can use other choices of W\ (H) and 
then update the resulting bond model using either a local 
algorithm [l9| or a worm-type algorithm 26]. 

Our second algorithm applies to the standard Eulcrian- 
subgraph model ([2]) with n > 1. Let us first observe, 
using the relation cycles = bonds — vertices + compo- 
nents, that we can replace n by n c ( H ^ [c(H] = cyclomatic 
number of H] in j2|) if we simultaneously replace Vij by 
%ij = Vij/n. Note also that the total number of faces 
(including the exterior face) in a bond configuration A 
equals its cyclomatic number plus 1. Therefore, we can 
attribute a weight n to each face. 

Algorithm 2: We again use two colors, but this time 
we color the faces of A (which are clusters of vertices of 
the dual graph G*V_In detail, the algorithm proceeds as 
follows (see again [25| for a proof of validity) : 

1) Independently for each face of A, color it a = 1 with 
probability l/n and a = 2 with probability 1 — 1/n. This 
decomposes the vertex set of G* as V\ U V% . 

2) On each edge of G* that does not lie entirely within 
V\ , place a bond. On each edge ij of G* that lies entirely 
within V\ and for which J e o~iO-j > 0, place a bond with 
probability 1 — e~ 2 l ,7,! L (Here e is the bond of G dual to 
ij, and J e is defined by x e — e~ 2l7c .) 

3,4) As in Algorithm 1. 

This algorithm is again easily extended to allowing 
more than one active color, in case n > 2, using the 
same rule as in Algorithm 1. 

Numerical results. We began by investigating the 
loop model ©/©on the honeycomb lattice, for which 
Nienhuis and Baxter 27 1 found the exact critical point 
when -2 < n < 2: rfnl = (2 + y/2 - n)" 1 / 2 . Nien- 
huis further observed [9|,|28( that the critical O(n) model 
with n > corresponds to a tricritical Potts model with 
q = n 2 . From Coulomb-gas theory, the leading and sub- 
leading thermal exponents and leading magnetic expo- 
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Vto 


Vti 


yho 




1.25 Exact 


1.83277 


0.88740 


1.93436 


1.38908 


Num. 


1.8327(1) 


0.887(2) 


1.9343(2) 1.3890(1) 


1.50 Exact 


1.78054 


0.74811 


1.91989 


1.40649 


Num. 


1.7805(1) 0.749(5) 


1.9198(1) 


1.4065(1) 


1.75 Exact 


1.70786 


0.55428 


1.90347 


1.43071 


Num. 


1.7080(1) 


0.55(1) 


1.9035(1) 


1.4307(1) 


2.00 Exact 


1.5 





1.875 


1.5 


Num. 


1.5000(1) 


-0.01(2) 1.8750(1) 


1.5001(1) 



Table I: A comparison of the critical 
by Monte Carlo simulations and those 



exponents determined 
predicted by (0 et seq. 



nent for the Potts model are known to be § 

3 5 -6 4<?-16 (0 + 2)0? + 6) 



yto 



<j 



yn 



ym 



8.9 



where g is the Coulomb-gas coupling and q = 
4cos 2 (7r<7/4); here g € [2,4] for critical and [4,6] for tri- 
critical. We also have a hull exponent [29( y H = (g+2) jg. 

We simulated the honeycomb-lattice loop model 
©/© at x = x c (n) for n = 1.25, 1.5, 1.75 and 2, us- 
ing Algorithm 1, for 14 linear lattice sizes L in the range 
4 < L < 1024. Periodic boundary conditions were ap- 
plied 0. We measured the observables N = number 
of occupied bonds, S' 2 = sum of squares of loop lengths, 
T>2 = sum of squares of face sizes, and M. = total Ising 
magnetization. From these we obtained the specific-heat- 
like quantity C = L~ d ({J\f 2 ) - (A/ - ) 2 ) and the Ising sus- 
ceptibility xis = L~ d (M 2 ) (here d = 2). 

Our numerical results for the static critical exponents 
are shown in Table [U We find empirically that xi s cx 
L 2 Vt0 -d^ L -d^ ^ const + i a «- d , C cx const + L 2 ^ 1 - d , 
L- d {V 2 ) cx L 2 ^- d and L- d {S' 2 ) cx L 2 y H ~ d . (We shall 
give elsewhere [2(| theoretical arguments for several of 
these.) The leading Potts exponents yto and yho are ab- 
sent in the spin and loop observables of the 0(n) model, 
but they can be seen in observables associated to the 
faces. To our knowledge, this is the first observation of 
yto, ym and yn in 0(n) models (see also 120]). 

We have also determined the integrated autocorrela- 
tion time r int for the observables Af, S 2 , T> 2 and A"! 2 . 
For 1 < n < 2, critical slowing-down is either entirely 
absent or very nearly absent. The Tint data for n = 2 
are shown in Fig. [T] and they strongly suggest that the 
dynamic critical exponent z is zero for n = 2. 

We next simulated model ([2|) on the square lattice, us- 
ing Algorithm 2, for lattices of linear size 4 < L < 512 
with periodic boundary conditions [3(J. The phase dia- 
gram of this model is largely unexplored, especially when 
n > 2. Phase transitions were located by analyzing the 
data for several Binder- type ratios. Two distinct phase 
transitions were found, which correspond to the ferro- 
magnetic and antiferromagnetic transitions in terms of 
the dual Ising spins (see Fig. [5]) . The three phases are, 
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1000 



Figure 1: Integrated autocorrelation times Tint for the 0(2) 
loop model on the honeycomb lattice versus lattice size L. 




Figure 2: Location of the ferro- and antiferromagnetic phase 
transitions of the standard Eulerian-subgraph model ((2} on 
the square lattice, as a function of n. 



respectively, dilute-dense, dense-dense and dense-dilute, 
where A-B means that the bond configuration is of type 
A and its complement is of type B (here "dense" = "per- 
colating"). The locations of the ferromagnetic critical 
point for n < 2 agree accurately with those given in 31 1 . 
Our data show that the lines of phase transitions extend 
to n > 2. We shall discuss elsewhere 25[ the nature of 
this new phase transition. The critical slowing-down is 
essentially absent for 1 < n < 2, but is strong for n > 2. 

We also applied Algorithm 1 to model (g]) on the square 
lattice. The results, along with a conjectured RG flow in 
the (n,v,t4) space, will appear elsewhere (25j |. 
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